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A  Monte  Carlo  ray  trace  radiation  model  is  presented  for  the  determination  of  radiative 
properties  of  A1203  particles  in  the  high  altitude  plume  of  a  solid  propellant  rocket.  A 
polydisperse  distribution  of  non-gray  particles  is  modeled  as  an  emitting,  absorbing  and 
scattering  medium  of  arbitrary  optical  thickness.  Strong  two-way  coupling  is  allowed 
between  radiation  and  flowfield  calculations,  where  the  gas  is  simulated  using  the  direct 
simulation  Monte  Carlo  method  and  particle  phase  properties  are  determined  using  a 
similar  Lagrangian  approach.  Effects  of  anisotropic  scattering  and  nozzle  searchlight 
emission  are  considered,  and  a  procedure  is  described  for  the  calculation  of  spectral 
radiance.  The  model  is  applied  to  the  simulation  and  radiation  analysis  of  the  freely 
expanding  plume  from  a  subscale  solid  rocket  motor,  and  various  flowfield  properties  are 
presented  and  discussed. 


I.  Introduction 

THE  analysis  and  prediction  of  radiation  signatures  from  solid  propellant  rocket  plumes  has  been  the  subject  of 
extensive  study  over  the  past  several  decades.1  Much  work  has  focused  on  properties  of  micron-scale  A1203 
particles  which  account  for  up  to  30%  of  the  mass  flow  through  the  nozzle,  and  which  tend  to  dominate  IR  radiative 
properties  within  the  plume.  Characteristics  of  the  particle  phase  are  extremely  complex,  and  are  generally  either 
poorly  understood  or  difficult  to  model  within  a  numerical  simulation.  Of  particular  interest  here  are  the  particle 
phase  properties  in  plume  flows  at  very  high  altitudes,  where  interactions  between  the  particles  and  gas  are  both 
important  and  not  well  suited  to  traditional  simulation  techniques. 

In  a  typical  high  altitude  exhaust  flow  from  a  solid  rocket  motor  (SRM),  liquid  aluminum  droplets  are  formed 
along  the  propellant  grain  surface,  and  undergo  a  complicated  process  of  combustion,  vaporization,  agglomeration, 
breakup,  and  crystallization  as  they  are  forced  by  the  expanding  gas  through  the  combustion  chamber  and  nozzle.2 
At  the  nozzle  exit,  the  particle  phase  consists  of  a  polydisperse  distribution  of  spherical  AFO,  particles,  which  range 
in  diameter  from  about  0.1  to  10  pm  depending  on  the  size  of  the  nozzle.  Larger  particles  exist  mainly  in  liquid  form 
within  the  nozzle  and  tend  to  develop  greater  temperature  and  velocity  lags  relative  to  the  gas.  The  smallest  particles 
are  rapidly  accelerated  and  cooled  by  the  surrounding  gas,  and  will  have  fully  solidified  into  some  combination  of 
metastable  gamma  phase  and  stable  alpha  phase  polycrystalline  structures  as  they  pass  through  the  nozzle  exit  plane. 

Within  the  plume  nearfield  region,  the  crystallization  process  becomes  increasingly  important  for  particles  of 
intermediate  size,  for  which  radiative  and  convective  cooling  are  balanced  by  a  heat  release  during  the  phase 
transition  process.  The  interaction  between  the  gas  and  particles  in  this  region  is  both  significant  and  complex,  as  the 
gas  develops  a  high  degree  of  thermal  nonequilibrium  during  the  rapid  expansion  downstream  of  the  nozzle  exit. 
Due  to  large  particle  phase  velocity  and  temperature  lags,  the  interphase  transfer  of  momentum  and  energy  may 
greatly  affect  particle  properties,  while  the  considerable  particle  mass  fraction  allows  bulk  gas  properties  to  be 
influenced  by  gas-particle  collisions.  Further  downstream,  interphase  momentum  and  energy  exchange  become 
negligible  as  the  gas  continues  to  expand.  The  largest  particles  may  begin  to  solidify  far  downstream  of  the  nozzle, 
where  radiative  heat  transfer  from  and  between  particles  tends  to  dominate  the  particle  phase  energy  balance. 
Depending  on  the  SRM  size  and  grain  composition,  the  farfield  plume  regions  will  likely  have  some  intermediate 
optical  thickness,  so  that  long-range  radiative  energy  exchange  within  the  particle  phase  may  significantly  influence 
particle  temperatures. 
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For  accurate  prediction  of  the  radiation  signature  from  a  high  altitude  SRM  plume  flow,  several  of  these 
processes  must  be  considered,  and  the  nonequilibrium  nature  of  the  gas  requires  a  simulation  approach  which  can 
overcome  limitations  of  traditional  computational  fluid  dynamics  (CFD)  techniques.  Various  complex  radiative 
properties  of  multiphase  particles  should  be  accounted  for,  including  a  temperature,  size,  and  wavelength 
dependence  in  spectral  emissivities.  Emission,  absorption,  and  scattering  may  all  be  important,  and  strong  two-way 
coupling  can  exist  between  the  radiation  field  and  particle  temperatures.  This  last  property  may  produce  large  errors 
in  the  traditional  post-processing  approach  to  plume  radiation  analysis.  With  these  goals  in  mind,  a  new  procedure  is 
proposed  for  the  simulation  and  radiation  analysis  of  high  altitude  SRM  plume  flows. 

Several  steps  in  this  procedure  have  been  extensively  documented  in  recent  papers,  so  are  only  discussed  here  in 
general  terms.  The  gas  is  simulated  using  the  Direct  Simulation  Monte  Carlo  (DSMC)  method  of  Bird,3  and  the  solid 
phase  is  modeled  using  a  similar  Lagrangian  approach  in  which  representative  particles  are  tracked  through  the 
computational  grid.  Momentum  and  energy  transfer  to  a  particle  from  the  surrounding  gas  is  modeled  through  a 
Green’s  function  technique  of  Gallis  et  al.,4  where  the  total  force  and  heat  transfer  rate  is  calculated  during  each  time 
step  by  summing  contributions  from  all  DSMC  gas  molecules  assigned  to  the  same  grid  cell.  The  reciprocal 
momentum  and  energy  transfer  from  the  particle  phase  to  the  gas  is  computed  using  a  probabilistic  approach,  in 
which  individual  interphase  collisions  are  modeled  as  involving  either  specular  reflection  or  diffuse  reflection  with 
full  thermal  accommodation.5  6  Several  additional  capabilities  and  physical  models  have  been  developed  and 
implemented  in  the  DSMC  code  MONACO.7  Mixed  discrete  and  continuous  particle  size  distributions  may  be  used, 
and  consideration  is  made  for  the  effects  of  particle  rotation,  nonspherical  particles,  the  breakup  of  liquid  droplets 
and  nonequilibrium  crystallization.  In  addition,  a  series  of  interphase  coupling  parameters8  are  utilized  to  increase 
computational  efficiency,  through  a  process  involving  the  automatic  determination  of  flowfield  regions  where 
calculations  for  momentum  and  energy  transfer  in  either  direction  may  be  avoided  with  negligible  impact  on  bulk 
flow  properties. 

In  the  following  sections,  a  particle  radiation  model  is  presented  for  use  with  the  above  procedures  in  the 
simulation  of  a  high  altitude  SRM  plume  flow.  First,  a  detailed  description  of  the  method  is  given,  including 
discussion  of  procedures  related  to  emission,  absorption,  scattering  and  two-way  coupling  between  radiation  and 
flowfield  characteristics.  In  addition,  techniques  are  outlined  for  consideration  of  nozzle  searchlight  emission  and 
the  calculation  of  spectral  radiance.  The  radiation  model  is  then  applied  to  the  plume  simulation  for  a  representative 
subscale  SRM  exhausting  into  a  vacuum.  Simulation  results  are  discussed,  and  sensitivities  of  radiation 
characteristics  to  various  input  parameters  are  evaluated.  In  particular,  we  consider  the  effects  of  nozzle  searchlight 
emission  and  anisotropic  scattering  on  spectral  radiance  and  the  net  radiative  energy  flux. 

II.  Radiation  Modeling  Procedure 

The  proposed  particle  radiation  model  uses  a  Monte  Carlo  Ray  Trace  (MCRT)  approach,  where  a  Lagrangian 
representation  is  used  to  track  large  groups  of  photons  through  the  computational  grid.910  As  a  first  step,  the  portion 
of  the  spectrum  of  interest  for  radiative  heat  transfer  -  wavelengths  of  roughly  0.5  to  5  pm  -  is  divided  into  a  series 
of  wave  number  bins.  Given  Nn  different  bins,  each  of  width  Aqt  and  centered  at  a  wave  number  q,,  a  large  number 
Nb  of  representative  “energy  bundles”  are  generated  once  every  several  time  steps  at  randomly  selected  source 
particles  throughout  the  grid.  Every  bundle  represents  some  quantity  of  radiative  power,  and  an  equal  number  of 
bundles  is  assigned  to  each  of  the  Nn  bins.  A  newly  generated  bundle  is  given  a  direction  of  propagation  according 
to  a  randomly  generated  unit  vector  u.  The  bundle  is  also  given  some  initial  power  Pb,  based  on  the  assigned  wave 
number  bin  and  the  properties  of  the  source  particle. 

Following  a  correlation  of  Reed  and  Calia1  based  on  Mie  theory  calculations  of  Plass,11  the  band-averaged 
spectral  emissivity  of  the  source  particle  at  bin  i  may  be  approximated  as 

=4k(Tp,r|i)Rpr|i  (1) 

where  Tp  and  Rp  are  the  temperature  and  radius  of  the  particle  respectively,  and  k  is  a  value  for  the  absorption  index 
of  AFOi  at  temperature  Tp  and  wave  number  q,.  Due  to  both  a  lack  of  experimental  data  and  an  extreme  sensitivity 
of  k  for  solid  particles  on  lattice  defects  and  impurities,1  we  neglect  any  dependence  of  k  on  the  particle  phase 
composition.  By  applying  Eq.  (1)  to  Planck’s  blackbody  function,  we  can  compute  the  source  particle  emissive 
power  Pp ;  within  the  wave  number  range  to  which  the  bundle  is  assigned: 

Pp,i  =  327T  cj;h  R3  Qj  (Tp )  (2) 
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where 


,  ktT.Wdl  (3) 

1  ■>  Arii  exp(hC()ri  /  kBT)  - 1 

The  symbol  c0  here  is  the  speed  of  light  through  a  vacuum,  h  is  Planck’s  constant,  and  kB  is  Boltzmann’s  constant. 
The  initial  power  of  the  bundle  is  then  determined  as 


N  N 

p  =^WnPn 
b  Nu  p  p 


(4) 


where  Np  is  the  total  number  of  representative  particles  in  the  grid  and  Wp  is  the  numerical  weight  of  the  source 
particle.  Values  of  Qj(T)  are  calculated  at  simulation  startup  for  each  wave  number  bin,  at  temperatures  T  for  which 
experimentally  determined  k(T,  ip)  values  are  available.  The  evaluation  of  Qj(Tp)  in  Eq.  (2)  is  then  performed 
through  linear  interpolation  to  the  particle  temperature. 

Once  created,  each  energy  bundle  is  tracked  through  the  grid  during  the  current  time  step  until  it  exits  through  an 
inflow,  outflow,  or  absorbing  wall  boundary.  As  an  energy  bundle  passes  through  a  cell  in  which  particles  are 
located  (or  have  been  located  during  previous  time  steps)  a  fraction  of  the  assigned  power  will  be  absorbed,  and 
there  is  some  probability  that  the  bundle  will  be  scattered.  Absorption  and  scattering  properties  in  the  cell  for  a  given 
wave  number  bin  are  determined  by  the  band-averaged  spectral  absorptance  a;  and  scattering  coefficient  a;. 
Consider  a  simulation  involving  Nspec  different  particle  species  j,  where  species  are  designated  according  to  the 
particle  radius  R,.  In  the  cell  of  interest,  each  species  has  an  average  temperature  T,  and  a  number  density  where 
averaging  is  performed  over  a  large  number  of  time  steps  for  both  Tj  and  nv  By  applying  Eq.  (1)  and  Kirchkoff  s  law 
to  the  definition  of  spectral  absorptance,  we  can  compute  a;  as 

Ns  ec 

«,  =  4mli  ^  k(Tj,rii )R3  n}  .  (5) 

j=i 


The  value  of  k(Tj,  rp)  for  each  species  is  found  by  interpolating  tabulated  k  values  to  the  temperature  Tj.  The 
corresponding  scattering  coefficient  C;  is  also  given  as  a  summation  over  all  particle  species: 

Ns  ec 

ai=7TE®i.jRi2«j  (6) 

j=l 


The  symbol  @y  in  Eq.  (6)  is  the  scattering  efficiency  factor  for  wave  number  bin  i  and  particle  species  j.  Values 
of  ©y  are  calculated  at  simulation  startup,  using  a  first-order  Mie  theory  approximation  of  Siegel  and  Howell.12 
Assuming  that  n;»k,  where  n,  is  the  real  part  of  the  index  of  refraction  for  the  particle  material  at  wave  number  q,, 
©y  may  be  given  as  a  function  of  the  nondimensional  parameter  xy=2jrr|iRj. 


4 

,i 3  xi. 


nr  -l 


n2  +2 


5 


V- 


nr  +2 


X; 


(7) 


While  Eq.  (7)  can  be  assumed  accurate  for  AI2O3  when  Xy«l,  it  greatly  over  predicts  ©y  for  larger  values  of  Xy.  To 
allow  for  use  with  a  wider  range  of  Xy  values,  we  impose  a  limiting  condition  ©y  <  2.  This  gives  relatively  good 
agreement  with  Mie  theory  calculations  of  Plass,11  and  avoids  the  detailed  calculations  required  to  find  an  exact  Mie 
theory  solution.  Note  that  we  neglect  here  any  dependence  of  n;  on  the  particle  temperature,  following  observations 
that  n,  values  for  AEO3  are  nearly  constant  over  a  wide  range  of  temperatures.11  While  experiments  have  shown 
some  increase  in  n,  with  particle  size,13  we  neglect  this  dependence  as  well  due  to  a  lack  of  available  experimental 
data. 
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When  an  energy  bundle  enters  a  given  cell,  the  total  distance  Dc  to  exit  the  cell  along  the  initial  trajectory  is 
determined  along  with  a;  and  a;  values  corresponding  to  the  assigned  wave  number  bin.  The  distance  Ds  to  a 
scattering  event  is  then  determined  by  evaluating  the  probability  Pns  that  the  bundle  will  not  have  been  scattered 
after  it  has  traveled  a  distance  Ds,  where 


dP^ 

dD 


=  —O:  P 


(8) 


Solving  for  Ds  and  setting  Pns  equal  to  a  random  number  9?e[0, 1 1.  we  find 

Ds=-— In  9?.  (9) 

If  DS>DC  then  the  particle  will  exit  the  cell  along  its  initial  trajectory.  Otherwise  the  bundle  will  be  scattered.  If 
scattered,  the  bundle  is  moved  a  distance  Ds  along  the  trajectory,  after  which  its  direction  is  reassigned  and  the 
procedure  is  repeated. 

The  anisotropic  nature  of  the  scattering  process  is  approximated  through  the  use  of  the  Henyey  Greenstein 
scattering  phase  function14 


0(0)  = 


1-g2 


47t(l  +  g2  -2gcos9) 


(10) 


The  free  parameter  g  in  Eq.  (10)  is  the  average  cosine  of  the  scattering  angle  0.  We  can  recover  the  corresponding 
distribution  function  /(9)  =  2ji<P(9)sin0  if  we  set 


jR  =  2 rcj  0(0  )sin0  dQ 


(11) 


for  a  random  number  9?  in  between  0  and  1.  Following  Eqs.  (10)  and  (11),  we  determine  0  through  the  formula 


cos  0  =  — 
2g 


g2+l- 


'  1-g2  A 

2g9t-g-l 


(12) 


The  unit  vector  u  in  the  final  direction  of  propagation  is  then  calculated  as 

u  =  u  cos  0  +  tj  sin  0  cos  ()>  + 1,  sin  0  sin  ()> 


(13) 


where  the  azimuthal  angle  (|>  is  assigned  a  random  value  in  [0,27t],  u  is  the  initial  direction,  and  the  unit  vectors  tj 
and  t2  are  given  by 


t,  = 


u  x  i 
u  x  i 


and  t2  =tj  xu  . 


For  convenience,  i  is  defined  here  as  the  unit  vector  along  the  x-coordinate  axis. 

The  procedure  involving  the  evaluation  of  Eqs.  (9),  (12)  and  (13)  is  repeated  until  the  bundle  exits  the  cell. 
Given  a  total  distance  Dt  which  the  bundle  has  traveled  through  the  cell,  the  assigned  power  Pb  is  then  reduced  by  a 
fraction  l-exp(-a.jDt)  to  account  for  the  effect  of  particle  phase  absorption  on  the  transmitted  radiation  intensity.9 

As  the  bundle  passes  through  the  cell,  the  power  AQabs  absorbed  by  an  individual  particle  of  radius  Rp  may  be 
given  as 
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(14) 


AQ 


abs 


^^Pb[l-exp(-a,Dt)] 

Ui  Vcell 


where  Vceu  is  the  cell  volume  and  Pb  is  the  initial  power  assigned  to  the  bundle  on  entering  the  cell.  It  follows  that 
the  contribution  of  the  bundle  to  the  direction-averaged  radiative  energy  flux  A q,  for  the  corresponding  wave  number 
bin  is 


A<?i  = 


— ^  [l -exp(-a;Dt)] 

Vcell 


(15) 


In  an  optically  thin  cell  where  a,Dt«l  the  evaluation  of  Eq.  (15)  may  result  in  a  large  subtractive  cancellation  error. 
To  correct  for  this  error,  when  a.iDt<10"5  we  calculate  Aq,  using  a  linearized  form  of  Eq.  (15) 

PhD, 

AO;=D^_L  (16) 

V 

v  cell 

based  on  a  Taylor  expansion  of  the  exponential  term.  Note  that  Eq.  (16)  gives  an  exact  solution  for  Aq,  in  regions 
outside  the  particle  domain  where  a;  is  zero.  During  each  time  step  for  which  energy  bundles  are  tracked  through  the 
grid,  energy  flux  values  in  every  cell  are  determined  by  summing  contributions  Aq,  from  all  bundles  which  pass 
through  the  cell.  The  resulting  values  are  then  averaged  over  a  large  number  of  time  steps  to  reduce  statistical 
scatter.  Once  the  flowfield  has  reached  steady  state  conditions,  averaging  may  be  performed  over  all  subsequent 
time  steps  during  which  radiation  calculations  are  made.  Otherwise  a  subrelaxation  technique  of  Sun  and  Boyd15  is 
used  for  the  time-averaging  procedure,  so  that  increased  weighting  may  be  applied  to  more  recent  time  steps. 

As  discussed  above,  strong  two-way  coupling  may  exist  between  flowfield  characteristics  and  plume  radiation. 
Radiative  heat  transfer  can  significantly  affect  particle  temperatures,  and  may  indirectly  influence  other  properties 
such  as  the  particle  phase  composition,  material  density,  and  the  rates  of  momentum  and  energy  transfer  between  the 
particles  and  gas.  To  account  for  the  effect  of  radiative  emission  and  absorption  on  particle  temperatures,  the 
temperature  Tp  of  every  representative  particle  is  modified  during  each  time  step  by 

ATp=—  Qrad  (17) 

mrcr 

where  At  is  the  time  step  interval,  mp  is  the  particle  mass,  cp  is  the  specific  heat  of  the  particle  material,  and  Qraij  is 
the  net  rate  of  radiative  heat  transfer  to  the  particle.  Note  that  the  particle  temperature  is  assumed  spatially  uniform, 
based  on  a  low  Biot  number  approximation  which  follows  from  the  small  particle  size  and  relatively  high  thermal 
conductivity  of  AfO>  The  radiative  heat  transfer  rate  is  then  calculated  as 

Nn  -. 

Qrad  =  4ttRp  ^  [ri^Tp ,  ri; )  -  8tt  Cq  h  Qj  (Tp )]  (18) 

i=i 

following  Kirchoff  s  law  and  Eqs.  (1)  and  (2).  The  above  symbol  q,  is  the  time-averaged  and  direction-averaged 
energy  flux  within  a  wave  number  bin  i,  for  the  cell  in  which  the  particle  is  located. 

Radiative  emission  from  within  the  nozzle  has  been  found  to  significantly  increase  radiance  intensity  near  the 
exit  plane.16  This  emission  is  generated  primarily  by  the  inner  nozzle  walls,  and  is  generally  termed  “searchlight 
emission”  under  past  assumptions  that  its  dominant  source  is  upstream  of  the  throat.  Depending  on  the  optical 
thickness  of  the  exhaust  flow,  emission  from  particles  within  the  nozzle  may  also  contribute  significantly.  As 
searchlight  emission  is  expected  to  influence  the  temperature  and  phase  composition  of  particles  in  the  plume,  a 
coupled  approach  to  radiation  and  flowfield  simulation  should  include  consideration  of  this  effect.  We  account  for 
searchlight  emission  through  the  generation  of  additional  energy  bundles  along  inflow  boundaries  at  the  nozzle  exit. 

Each  inflow  boundary  on  the  exit  plane  is  represented  as  a  blackbody  wall  at  some  characteristic  temperature  Tw. 
Along  every  cell  face  located  on  an  inflow  boundary,  N,  new  bundles  are  generated  during  each  time  step  for  which 
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radiation  calculations  are  performed.  Each  bundle  is  randomly  assigned  to  a  wave  number  bin  i,  and  is  given  an 
initial  direction  u  such  that 


unf  =  cos  0  =  SJ?12 


where  nf  is  the  inward  normal  unit  vector  at  the  cell  face,  0  is  the  zenithal  angle  of  u  relative  to  the  face,  and  5?  is  a 
random  number  in  [0,1].  The  initial  power  Pb  is  then  determined  from  Eq.  (19)  as  the  product  of  the  face  area  Af,  a 
weighting  factor  Nn/Nr,  and  the  integral  of  Planck’s  function  over  the  corresponding  wave  number  bin. 


Pb  =  2nc20  h  Af 


N 


r|3dr| 


Nf  Jati.  exp(hc0r| / kBTw )  - 1 


(19) 


For  calculations  of  plume  radiance,  one  or  more  simulated  sensors  are  positioned  somewhere  outside  the  grid 
domain.  Consider  a  sensor  with  surface  area  As,  outward  unit  normal  ns,  and  angular  resolution  defined  by  the 
zenithal  angle  to.  When  each  energy  bundle  exits  the  grid,  we  determine  whether  it  will  intersect  the  sensor  surface 
along  a  trajectory  given  by  the  unit  vector  u  for  the  bundle  direction.  If  an  intersection  occurs  and  the  condition 


-u  ns  >  cos  to 


is  met,  then  the  power  Pb  assigned  to  the  bundle  is  added  to  the  total  absorbed  power  IP,  for  the  corresponding  wave 
number  bin  during  the  current  time  step.  An  instantaneous  value  of  the  band-averaged  spectral  radiance  I;  may  then 
be  calculated  as  the  ratio  of  ZP;  to  the  product  of  the  sensor  area,  solid  angle  for  absorbed  radiation,  and  wavelength 
range  of  the  bin.  We  find 


I-  =- 


IP 


2tc(1  -coseo)As  AXt 


(20) 


where 


AX;  =  ' 


Atf 


nr-W 


Values  of  I;  for  each  bin  are  averaged  over  a  large  number  of  time  steps  to  reduce  scatter,  where  sampling  is 
begun  only  after  the  flowfield  has  reached  steady  state  conditions.  Note  that  the  rate  of  statistical  convergence  in  the 
resulting  radiance  values  will  be  roughly  equal  for  every  bin.  While  a  much  faster  convergence  rate  is  possible  if  we 
apply  a  reverse  Monte  Carlo  method  involving  the  generation  of  additional  energy  bundles  at  the  sensor,14  the 
procedure  proposed  here  allows  for  relatively  fast  convergence  without  adding  much  to  the  complexity  of  the 
radiation  model  implementation. 


III.  Plume  Flow  Test  Case 

Following  implementation  in  the  DSMC  code  MONACO,7  the  procedures  described  above  are  applied  to  an 
axisymmetric  test  case  for  the  plume  flow  from  a  subscale  SRM  expelling  into  a  vacuum.  Inflow  and  boundary 
conditions  are  identical  to  those  described  in  a  previous  paper,8  as  are  all  physical  models  used  in  the  flowfield 
simulation.  We  use  a  rectangular  grid  geometry  extending  from  the  nozzle  exit  plane  100  m  downstream  and  40  m 
radially  outward.  The  nozzle  exit  diameter  is  7.85  cm,  and  exit  plane  data  for  both  the  particles  and  gas  are  taken 
from  Anfimov  et  al.17  based  on  simulated  nozzle  flow  characteristics  for  a  Star-27  motor  with  30%  particle  mass 
loading.  The  gas  is  a  mixture  of  N2,  H2  and  CO,  and  the  variable  hard  sphere  (VHS)  collision  model3  is  used  for 
collisions  within  the  gas  phase.  At  the  exit  plane,  the  gas  is  assigned  a  bulk  speed  of  3113  m/s,  a  temperature  of 
1433  K  and  a  density  of  0.01 1  kg/m3,  with  mole  fractions  of  0.38  for  H2  and  0.3 1  for  both  N2  and  CO. 

The  particle  phase  has  a  discrete  size  distribution,  with  seven  different  species  ranging  in  diameter  from  0.3  to  6 
pm.  Due  to  a  lack  of  available  flowfield  information  at  the  nozzle  exit,  particle  properties  given  by  Anfimov  et  al.17 
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are  assumed  uniform  over  the  exit  plane.  A  rough  estimate  of  the  initial  solid  mass  fraction  for  each  particle  size  is 
found  by  taking  the  difference  between  the  nucleation  temperature  (1930  K)  for  homogeneous  crystallization  of 
AI2O3  and  the  assigned  particle  temperature,  then  multiplying  this  by  the  ratio  of  the  specific  heat  to  the  latent  heat 
of  fusion  for  liquid  A1203.  Particle  properties  at  the  nozzle  exit  are  given  in  Table  1. 

The  thermal  accommodation  coefficient  on  the  particle  surface  is  set  to  0.9,  so  that  90%  of  inteiphase  collisions 
involve  diffuse  reflection  with  full  accommodation  to  the  particle  temperature,  while  the  remaining  10%  involve 
specular  reflection.  Interphase  momentum  and  energy  transfer  is  computed  using  the  two-way  coupling  approach 
described  above,  and  the  crystallization  of  liquid  A1203  droplets  is  considered  using  a  model  for  nonequilibrium 
phase  change.8  The  phase  change  model  accounts  for  a  temperature  dependence  in  the  crystallization  rate  and  the 
associated  heat  release,  while  neglecting  the  gamma-to-alpha  transition  for  solid  A1203  and  the  density  variation 
between  different  phases. 

Values  of  the  particle  phase  absorption  index  k  are  taken  from  experimental  data  of  Konopka,  Reed  and  Calia.18 
Based  on  this  data,  we  use  10  wave  number  bins  corresponding  to  the  mid-IR  wavelength  range  from  1.3  to  4.5  pm. 
Of  the  two  SRM  exhaust  flows  from  which  particles  were  collected  and  investigated  by  these  authors,  the  second 
(rocket  2)  was  found  to  give  k  values  more  in  line  with  other  experimental  data  and  correlations  in  the  literature.19 
Values  calculated  from  SEM  measurements  for  the  second  flow  are  therefore  used  here.  The  real  part  n,  of  the 
particle  index  of  refraction  is  computed  as  a  function  of  r),  through  a  correlation  given  by  Duval  et  al.,1 


:  (0.9904  +  2.02  xl(T5Tp) 


1  +  - 


1.024 


1.058 


5.281 


-1I/2 


1-0.00376TV  l-0.01225ri?  1-321.4t|? 


(21) 


where  the  weak  temperature  dependence  is  neglected  by  assuming  a  particle  temperature  Tp  of  2000  K. 

Following  Reed  et  al.,16  we  set  the  average  cosine  of  the  scattering  angle  to  g  =  0.5  and  use  an  effective 
temperature  of  Tw  =  1300  K  for  searchlight  emission  at  the  nozzle  exit.  The  unstructured  grid  consists  of 
approximately  40,000  triangular  cells  scaled  roughly  according  to  the  local  mean  free  path.  So  that  calculations  may 
be  performed  in  a  reasonable  amount  of  time  on  available  computing  resources,  cells  are  about  10  to  50  times  larger 
than  the  mean  free  path  in  a  small  region  just  beyond  the  nozzle  exit.  While  this  is  a  violation  of  the  standard  rule  of 
thumb  for  cell  size  in  DSMC,  we  can  tolerate  the  associated  errors  with  an  understanding  that  the  simulation 
performed  here  is  intended  less  to  represent  one  particular  flow  with  maximum  accuracy  than  to  show  the  general 
radiation  characteristics  of  a  typical  small-scale  SRM  plume  flow  at  high  altitude. 

About  two  million  DSMC  gas  molecules  and  200,000  representative  solid  particles  are  tracked  through  the  grid 
during  a  typical  time  step  at  steady  state.  Sampling  is  performed  once  every  five  time  steps,  when  1600  energy 
bundles  are  generated  and  moved  through  the  grid.  To  reach  the  desired  level  of  statistical  scatter  in  sampled 
radiation  and  bulk  flow  properties,  the  simulation  is  run  for  100  hours  on  four  processors  in  a  2.4  GHz  AMD  Athlon 
cluster.  Calculations  are  divided  roughly  evenly  across  all  processors  through  domain  decomposition,  and  buffer 
arrays  are  used  to  allow  energy  bundles  to  rapidly  move  across  multiple  task  domains  while  restricting  information 
exchange  between  neighboring  tasks.  Selected  simulation  results  are  shown  in  Figs.  1  through  8. 

IV.  Simulation  Results 

Mass  density  contours  for  both  particles  and  gas  are  shown  in  Fig.  1.  Note  first  that  particles  are  only  found  in 
roughly  half  of  the  simulation  domain,  as  the  maximum  divergence  angle  for  the  particle  phase  is  restricted  by 
particle  mass.  All  particles  are  given  initial  trajectories  between  0  and  15  degrees  off  the  centerline  at  the  nozzle  exit 
plane,  where  the  trajectory  angle  for  each  particle  scales  linearly  with  distance  from  the  axis.  In  the  plume  nearfield 
region  just  beyond  the  nozzle  exit,  particles  are  forced  outward  from  the  centerline  by  the  expanding  gas.  The  radial 
acceleration  of  an  individual  particle  in  this  region  will  vary  as  the  inverse  of  the  particle  diameter,  so  that  in  a  given 
radial  plane  smaller  particles  will  be  found  over  a  range  which  extends  further  from  the  axis.  As  shown  in  the  top 
half  of  the  figure,  this  results  in  a  gradual  decrease  in  particle  mass  density  with  distance  from  the  axis  due  to  the 
presence  of  a  range  of  particle  sizes. 

As  the  dimensions  of  the  simulation  domain  are  several  orders  of  magnitude  greater  than  the  nozzle  exit  radius, 
particle  characteristics  observed  in  Fig.  1  reflect  only  trends  in  the  farfield  region,  where  momentum  and  energy 
transfer  between  the  particles  and  gas  can  be  assumed  negligible.  The  particles  move  along  nearly  straight 
trajectories  far  from  the  nozzle,  so  the  contour  lines  shown  projecting  from  the  nozzle  exit  are  straight  if  we  neglect 
effects  of  statistical  scatter  and  interpolation  between  cell  centers.  The  continuous  reduction  in  particle  mass  density 
with  downstream  distance  is  due  to  the  divergence  of  particle  trajectories,  and  the  restricted  maximum  divergence 
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angle  is  shown  in  the  figure  to  result  in  centerline  values  of  the  particle  phase  mass  fraction  (the  ratio  of  particle 
mass  density  to  total  mass  density  for  both  particles  and  gas)  which  increase  with  distance  from  the  nozzle.  At  the 
point  on  the  centerline  100  m  downstream  of  the  nozzle  exit,  the  particle  mass  fraction  is  about  2.2  times  the  initial 
value  of  0.3. 

Figure  2  shows  contours  of  average  temperature  for  particles  of  diameter  0.4  pm  and  4  pm.  The  temperature  of 
0.4  pm  particles  is  generally  observed  to  decrease  with  downstream  distance  but  increase  with  distance  from  the 
axis.  The  first  trend  may  be  attributed  to  the  dominance  of  radiative  heat  loss  in  the  particle  energy  balance  through 
nearly  the  entire  simulation  domain,  while  the  second  trend  is  primarily  a  result  of  the  radial  variation  in  gas  density 
within  the  plume  nearfield  region.  As  gas  streamlines  diverge  in  this  region  just  downstream  of  the  nozzle  exit,  the 
density  of  the  gas  will  decrease  more  rapidly  along  streamlines  further  from  the  centerline.  A  lower  gas  density 
corresponds  to  a  reduction  in  the  rate  of  convective  heat  transfer  from  particles  to  the  surrounding  gas,  so  we  can 
expect  a  smaller  temperature  drop  in  the  plume  nearfield  for  particles  which  pass  through  the  nozzle  exit  plane 
further  from  the  centerline.  This  is  likely  the  dominant  mechanism  for  the  increase  in  particle  temperature  with 
distance  from  the  axis. 

Another  possible  contributor  to  this  trend  is  the  fact  that  particles  further  from  the  axis  may  absorb  more 
radiative  energy  from  nozzle  searchlight  emission.  The  plume  optical  thickness  scales  with  particle  phase  mass 
density,  and  this  density  is  shown  in  Fig.  1  to  decrease  with  distance  from  the  axis,  so  a  corresponding  increase  in 
the  magnitude  of  long  range  radiative  intensity  from  nozzle  emission  should  result  in  some  increase  in  absorption 
among  particles  far  from  the  centerline.  Flowever,  this  trend  is  found  to  have  only  a  very  small  effect  on  particle 
temperatures,  and  is  countered  by  an  increase  in  short  range  radiative  heat  transfer  between  particles  near  the  axis. 
We  expect  short  range  radiative  transfer  to  be  the  dominant  mechanism  by  which  radiation  causes  a  radial  variation 
in  particle  temperatures,  so  that  the  net  effect  of  coupling  the  radiation  model  to  the  flowfield  simulation  should  be  a 
slight  reduction  in  the  temperature  increase  for  0.4  pm  particles  with  distance  from  the  axis.  Note  that  a  significant 
temperature  drop  is  observed  for  these  particles  at  points  very  close  to  the  axis.  This  is  likely  the  result  of  an 
unphysical  decrease  in  convective  heat  transfer  between  the  particles  and  gas,  which  follows  from  a  reduction  in  the 
average  number  of  DSMC  gas  molecules  in  cells  along  the  centerline  just  downstream  of  the  nozzle  exit. 

Similar  trends  as  described  above  are  found  within  50  m  of  the  nozzle  exit  for  4  pm  diameter  particles,  as  shown 
in  the  lower  half  of  Fig.  2.  Average  temperatures  for  these  particles  in  the  first  50  m  downstream  of  the  nozzle  are 
several  hundred  degrees  higher  than  for  the  0.4  pm  particles.  This  may  be  explained  by  the  fact  that  the  drop  in 
particle  temperature  within  the  nozzle  and  the  nearfield  plume  regions  scales  roughly  with  the  inverse  of  the  particle 
diameter,  following  a  balance  between  particle  heat  capacity  and  the  heat  transfer  rate.  As  all  particles  enter  the 
nozzle  at  similar  temperatures,  we  can  generally  expect  an  increase  in  temperature  with  particle  size  throughout  the 
plume. 

Radiative  heat  transfer  is  shown  in  the  figure  to  result  in  a  continuous  downstream  reduction  in  temperature  for  4 
pm  particles,  up  to  a  narrow  region  about  50  m  away  from  the  nozzle  where  these  particles  reach  a  temperature  of 
1930  K.  In  our  model  for  nonequilibrium  phase  change,  this  is  specified  as  the  nucleation  temperature  at  which  a 
crystallization  front  forms  uniformly  over  the  surface  of  an  initially  liquid  particle.  Once  formed,  this  front 
progresses  toward  the  particle  center  at  a  temperature  dependent  rate,  so  long  as  the  particle  temperature  remains 
below  the  equilibrium  melting  temperature  of  2327  K.  As  the  heat  of  formation  for  liquid  AFO^  is  released  during 
crystallization,  particles  will  experience  a  rapid  temperature  increase  at  the  initiation  of  the  phase  change  process. 
The  rate  of  this  increase  will  be  reduced  further  downstream  as  the  particle  temperature  rises,  so  the  crystallization 
front  progresses  more  slowly,  and  as  the  front  area  decreases,  so  a  smaller  volumetric  phase  change  rate  will  exist 
for  a  given  velocity  of  progression  toward  the  particle  center.  Ultimately,  as  the  particle  becomes  completely 
solidified,  we  expect  the  particle  temperature  to  again  begin  to  decrease  with  downstream  distance  due  to  radiative 
heat  loss.  For  4  pm  particles,  the  temperature  reduction  following  crystallization  should  occur  beyond  the 
downstream  limit  of  the  simulation  domain,  so  this  trend  is  not  shown  in  Fig.  2.  Note  that  the  phase  change  process 
for  0.4  pm  particles  is  completed  primarily  within  the  nozzle,  so  no  similar  temperature  jump  is  shown  in  the  upper 
half  of  the  figure. 

In  Fig.  3,  particle  temperature  profiles  are  shown  along  the  radial  plane  100  m  downstream  of  the  nozzle  exit. 
Average  temperatures  for  particles  of  four  different  sizes  are  plotted  as  a  function  of  distance  from  the  axis,  and  for 
comparison,  corresponding  temperature  profiles  are  shown  for  a  simulation  in  which  the  radiation  model  is  disabled. 
First  consider  the  temperature  profiles  for  the  latter  case:  For  all  particle  sizes  shown,  temperatures  are  found  to 
increase  with  distance  from  the  axis  due  to  the  radial  variation  in  gas  density  within  the  plume  nearfield,  as 
discussed  above.  The  temperature  is  also  found  to  generally  increase  with  particle  size,  so  that  the  lowest 
temperatures  are  observed  for  the  smallest  particles  shown,  while  the  highest  temperatures  occur  for  the  largest 
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particles.  An  important  exception  to  this  trend  is  observed  for  0.6  pm  diameter  particles,  which  experience 
temperatures  significantly  higher  than  those  of  the  much  larger  4  pm  particles.  This  can  be  explained  by 
characteristics  of  the  crystallization  process.  While  4  pm  particles  never  undergo  sufficient  convective  cooling  to 
reach  the  nucleation  temperature  of  1930  K  at  which  crystallization  may  begin,  the  0.4  pm  and  0.6  pm  particles 
pass  through  the  nozzle  exit  plane  below  the  melting  temperature  and  in  a  state  of  partial  solidification.  The  entire 
phase  transition  process  for  0.4  pm  particles  takes  place  within  and  just  beyond  the  nozzle,  where  convective  heat 
loss  dominates  the  particle  energy  balance  and  allows  for  temperatures  several  hundred  degrees  below  the  melting 
temperature  of  2327  K.  In  contrast,  the  phase  change  process  for  0.6  pm  particles  occurs  at  a  slower  rate,  and 
continues  well  beyond  the  nearfield  plume  region  where  the  associated  heat  release  may  be  balanced  by  convective 
heat  transfer  to  the  gas.  As  a  result,  the  temperature  of  these  particles  approaches  a  farfield  limit  which  is 
significantly  greater  than  that  of  both  larger  and  smaller  particles. 

In  comparing  temperature  profiles  between  the  cases  with  and  without  activation  of  the  radiation  model,  we  find 
that  radiative  heat  loss  accounts  for  a  roughly  120  K  decrease  in  temperatures  for  0.4  pm  particles  and  a  200  K 
decrease  for  0.6  pm  particles.  Among  4  pm  particles  however,  radiative  heat  transfer  is  found  to  increase 
temperatures  100  m  downstream  of  the  nozzle  by  about  240  K.  As  above,  this  trend  may  be  explained  by  the  phase 
change  process:  Convective  cooling  within  the  nozzle  and  plume  nearfield  regions  is  not  sufficient  to  reduce  the 
temperature  of  4  pm  particles  to  below  the  nucleation  temperature,  but  the  addition  of  radiative  cooling  results  in 
farfield  particle  temperatures  at  which  crystallization  fronts  may  form.  Because  the  rate  of  heat  addition  for  a 
particle  undergoing  crystallization  is  far  greater  than  the  rate  of  radiative  heat  loss,  these  particles  will  experience  a 
rapid  temperature  increase  as  shown  in  Fig.  2. 

For  6  pm  particles,  a  comparison  of  temperature  profiles  shows  an  even  more  complicated  trend.  Radiation  is 
found  to  uniformly  reduce  particle  temperatures,  but  the  magnitude  of  this  reduction  varies  greatly  with  distance 
from  the  axis.  This  trend  may  be  explained  as  follows:  Radial  variation  in  nearfield  convective  heat  transfer  results 
in  an  increase  in  temperature  for  these  particles  with  distance  from  the  axis,  as  described  above.  This  means  that 
particles  further  from  the  axis  must  experience  radiative  heat  loss  over  a  longer  time  period  before  their  temperatures 
have  been  sufficiently  reduced  for  crystallization  to  begin.  The  temperature  rise  associated  with  crystallization 
therefore  occurs  further  downstream  for  particles  at  a  greater  distance  from  the  axis.  At  the  radial  plane  100  m 
downstream  of  the  axis,  6  pm  particles  close  to  the  centerline  have  begun  the  phase  change  process,  while  the 
fraction  of  particles  on  which  crystallization  fronts  have  formed  is  generally  found  to  decrease  with  radial  distance. 
This  reduction  in  liquid  mass  fraction  with  radial  distance  corresponds  to  the  temperature  drop  observed  in  Fig.  3  for 
6  pm  particles  far  from  the  axis. 

Figure  4  shows  the  centerline  variation  in  the  magnitude  of  mean  convective  and  radiative  heat  transfer  rates,  as 
calculated  per  particle  and  averaged  over  all  particle  sizes.  Both  rates  are  negative  through  the  entire  simulation 
domain,  so  that  particles  throughout  the  plume  are  losing  energy  to  their  surroundings  through  both  convective  and 
radiative  heat  transfer.  The  radiative  heat  transfer  rate  is  shown  to  be  relatively  constant  with  downstream  distance, 
due  to  an  overall  gradual  variation  in  particle  temperatures  and  the  fact  that  radiative  absorption  is  found  to  have  a 
comparatively  small  effect.  Most  of  the  variation  in  radiative  transfer  observed  on  the  plot  is  the  result  of  statistical 
scatter,  due  to  the  small  number  of  representative  particles  which  pass  through  cells  bordering  the  axis.  However,  a 
significant  increase  in  the  radiative  transfer  rate  is  found  about  50  m  downstream  of  the  nozzle  exit.  This  may  be 
attributed  to  a  temperature  jump  associated  with  the  onset  of  phase  change  in  4  pm  particles,  which  account  for 
about  60%  of  the  total  AFOj  mass  within  the  plume. 

In  contrast  to  the  relatively  uniform  rate  of  radiative  heat  transfer,  the  mean  convective  heat  transfer  rate  is 
shown  to  decrease  rapidly  with  downstream  distance.  The  spatial  variation  in  convective  heat  transfer  is  found  to 
occur  mainly  as  a  result  of  a  downstream  decrease  in  gas  density.  As  the  gas  approaches  a  free  molecular  state  far 
from  the  nozzle,  the  gas  density  will  decrease  along  the  centerline  as  the  inverse  square  of  the  distance  from  some 
point  near  the  nozzle  exit.  The  farfield  convective  heat  transfer  rate  will  therefore  have  the  same  inverse  square  of 
distance  variation,  as  is  shown  in  Fig.  4.  A  convenient  definition  of  the  plume  nearfield  region  in  a  freely  expanding 
SRM  plume  flow  is  the  range  beyond  the  nozzle  exit  where  the  particle  energy  balance  is  dominated  by  convective 
heat  transfer  to  the  gas.  By  this  definition,  the  nearfield  region  extends  along  the  centerline  about  1  m  (or  13  nozzle 
exit  diameters)  downstream  of  the  nozzle  exit,  beyond  which  radiative  emission  becomes  the  dominant  mechanism 
for  energy  transfer  between  a  particle  and  its  surroundings. 

In  Fig.  5,  contours  are  shown  for  the  net  direction-averaged  radiative  energy  flux.  The  energy  flux  is  greatest  at 
the  nozzle  exit,  due  to  the  corresponding  maximum  in  particle  mass  density  and  a  reduction  in  intensity  of 
searchlight  emission  with  distance  from  the  nozzle.  Particles  are  modeled  through  Eq.  (1)  as  volumetric  emitters,1  so 
the  magnitude  of  radiative  energy  flux  should  scale  roughly  with  the  local  particle  mass  density.  As  this  density 
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decreases  with  downstream  distance  due  to  the  divergence  of  particle  trajectories,  a  continuous  reduction  in  radiative 
energy  flux  will  occur  in  the  axial  direction. 

A  reduction  in  energy  flux  is  also  observed  in  the  radial  direction,  particularly  outside  the  region  where  particles 
are  found.  This  follows  primarily  from  the  fact  that  in  an  axisymmetric  simulation  the  ratio  of  cell  volume  to  the 
projected  area  on  the  grid  will  scale  with  distance  from  the  axis.  The  average  number  of  energy  bundles  which  pass 
through  a  cell  is  proportional  to  the  projected  area  of  that  cell,  and  the  contribution  of  each  bundle  to  the  energy  flux 
is  proportional  to  the  inverse  of  the  cell  volume.  We  therefore  expect  the  net  energy  flux  through  cells  outside  the 
particle  domain  to  scale  approximately  with  the  inverse  of  the  distance  from  the  axis,  as  is  shown  in  the  figure. 

Figure  6  shows  a  comparison  of  net  radiative  energy  flux  contours  between  simulations  with  and  without 
activation  of  the  model  for  nozzle  searchlight  emission.  Due  to  the  point-source  nature  of  this  effect,  searchlight 
emission  is  found  to  have  relatively  little  influence  on  energy  flux  magnitudes  far  downstream  of  the  nozzle.  For 
clarity  we  therefore  restrict  the  figure  to  an  axial  range  of  10  m  from  the  nozzle  exit.  The  presence  of  searchlight 
emission  results  in  a  100%  increase  in  radiative  energy  flux  at  the  point  along  the  centerline  0.5  m  from  the  nozzle 
exit.  At  points  1,  2,  5  and  10  m  downstream,  we  find  corresponding  increases  of  about  45%,  13%,  11%  and  9%, 
respectively.  An  equivalent  contour  plot  to  evaluate  the  impact  of  anisotropic  scattering  is  included  as  Fig.  7.  The 
top  half  of  this  figure  shows  contours  of  net  direction-averaged  energy  flux  for  the  base  simulation,  as  described 
above,  where  we  apply  the  Flenyey  Greenstein  scattering  phase  function.  The  lower  half  of  the  figure  is  taken  from  a 
simulation  where  isotropic  scattering  is  used,  so  that  each  scattered  energy  bundle  is  randomly  assigned  a  new 
direction  of  propagation  with  no  dependence  on  the  initial  direction.  We  find  little  if  any  measurable  effect  of  the 
scattering  model  on  energy  flux  contours,  as  the  small  differences  observed  between  the  upper  and  lower  halves  of 
Fig.  7  are  primarily  a  result  of  statistical  fluxuations  due  to  the  probabilistic  nature  of  the  radiation  model. 

Spectral  radiance  is  calculated  at  a  simulated  sensor  centered  at  a  point  20  cm  upstream  and  4  cm  radially 
outward  from  the  intersection  of  the  central  axis  with  the  nozzle  exit  plane.  The  sensor  has  a  circular  face  of  area  50 
cm2  and  a  surface  normal  vector  inclined  4°  from  the  axis.  The  angular  resolution  is  4°,  so  that  the  boundary  of  the 
conical  viewing  area  intersects  the  grid  plane  along  a  line  which  is  parallel  to  the  axis.  Values  of  spectral  radiance 
are  plotted  in  Fig.  8,  where  the  solid  line  denotes  the  base  simulation  and  dashed  lines  correspond  to  simulations  for 
which  either  searchlight  emission  or  the  anisotropic  scattering  model  is  disabled.  All  three  lines  generally  follow  the 
gray  body-like  trend  described  by  Reed  and  Calia.1 

As  expected,  we  find  that  searchlight  emission  produces  an  increase  in  radiance  values  upstream  of  the  nozzle. 
Flowever,  the  magnitude  of  this  increase  is  very  small  compared  to  the  increase  in  radiative  energy  flux  near  the 
nozzle  shown  in  Fig.  6.  As  discussed  by  Reed  et  al., 16  the  relatively  small  dependence  of  measured  radiance  on 
searchlight  emission  may  be  attributed  to  a  strong  preference  for  forward  scattering  among  A1203  particles.  Any 
nozzle  emission  absorbed  by  the  sensor  must  be  scattered  off  one  or  more  particles  at  a  large  angle,  for  which  the 
scattering  efficiency  is  very  low.  Searchlight  emission  will  therefore  have  a  far  smaller  effect  on  radiance  upstream 
of  the  nozzle  than  on  radiance  measured  at  a  point  downstream  of  the  exit  plane.  If  the  angular  dependence  in  the 
scattering  model  is  disabled,  then  we  expect  to  find  a  significant  increase  in  radiance  upstream  of  the  nozzle  due  to 
searchlight  emission.  This  trend  is  shown  in  Fig.  8,  where  a  jump  in  spectral  radiance  of  up  to  20%  is  found  when 
scattering  is  assumed  to  be  isotropic. 


V.  Conclusions 

A  radiation  model  was  presented  for  the  AI2O3  particle  phase  in  SRM  plume  flows  at  high  altitude.  The  model 
includes  capabilities  for  strong  two-way  coupling  between  radiation  and  flowfield  calculations,  and  accounts  for  the 
influence  of  emission  and  absorption  on  particle  temperatures.  No  limitations  are  imposed  on  optical  thickness  of  the 
flowfield,  and  band-averaged  particle  radiative  properties  are  allowed  to  vary  as  an  arbitrary  function  of  both 
temperature  and  wavelength.  Procedures  were  described  for  the  calculation  of  spectral  radiance  and  consideration  of 
effects  associated  with  nozzle  searchlight  emission.  Following  application  of  the  model  in  the  axisymmetric 
simulation  of  a  subscale  plume  flow,  several  radiation  and  flowfield  characteristics  were  evaluated.  A  complex 
relationship  was  shown  between  particle  size  and  temperature  in  farfield  plume  regions,  due  in  part  to  the  interaction 
between  radiative  heat  loss  and  crystallization  of  liquid  AFO3.  Searchlight  emission  was  found  to  significantly  affect 
the  radiative  energy  flux  through  a  large  region  beyond  the  nozzle  exit,  while  consideration  of  the  angular 
dependence  in  particle  scattering  results  in  an  intensity  reduction  upstream  of  the  nozzle. 

Due  to  a  lack  of  available  experimental  data  for  comparison,  the  overall  accuracy  of  our  radiation  model  in  the 
simulation  of  a  high  altitude  SRM  plume  flow  could  not  be  quantitatively  evaluated.  Flowever,  we  expect  that  the 
single  greatest  error  source  is  the  selection  of  appropriate  values  for  the  particle  absorption  index,  particularly  at  low 
temperatures  where  much  of  the  particle  material  has  solidified.  As  discussed  by  Reed  and  Calia,1  the  absorption 
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index  of  solid  AI2O3  particles  is  primarily  an  extrinsic  property,  and  is  a  strong  function  of  the  concentration  of 
impurities  within  the  lattice  structure.  As  impurity  concentrations  will  vary  greatly  between  different  SRM  exhaust 
flows,  and  between  different  regions  in  the  same  flow,  the  identification  of  reasonably  accurate  values  for  the 
absorption  index  of  solid  A1203  becomes  extremely  difficult.  Based  on  experimental  measurements18  we  can  expect 
these  values  to  be  accurate  to  about  one  order  of  magnitude,  so  that  radiative  heat  transfer  rates  and  energy  flux 
values  far  from  the  nozzle  may  vary  on  this  same  scale. 

Another  potential  source  of  significant  error  is  the  difficulty  in  accurately  determining  particle  properties  at  the 
nozzle  exit.  As  a  result  of  complicated  physical  processes  which  occur  within  the  combustion  chamber  and  nozzle, 
such  as  particle  agglomeration,  breakup,  combustion,  and  turbulent  dispersion,  the  accurate  numerical  calculation  of 
particle  properties  at  the  exit  plane  may  be  very  difficult.  Other  difficulties  are  associated  with  the  experimental 
measurement  of  particle  phase  characteristics  along  the  nozzle  exit,  and  no  sufficiently  detailed  experimental  data 
could  be  found  in  the  open  literature.  We  must  therefore  rely  on  simplified  nozzle  flow  simulations  which  neglect 
many  of  the  physical  phenomena  expected  in  such  flows,  so  that  a  loss  of  accuracy  will  result  in  the  calculation  of 
particle  properties  within  the  plume.  Other  possible  sources  of  significant  error  in  the  radiation  model  include  the 
lack  of  consideration  for  IR  absorption  and  emission  involving  exhaust  gas  species,  approximations  used  to  compute 
particle  emissivity  and  scattering  coefficients,  and  an  unphysical  reduction  in  farfield  energy  fluxes  and  radiative 
heat  transfer  due  to  the  truncation  of  the  simulation  domain. 

While  the  radiation  model  presented  here  was  developed  for  application  to  high  altitude  SRM  plume  flows,  the 
above  procedures  may  be  applied  to  any  simulation  involving  a  Lagrangian  representation  of  micron-scale  A1203 
particles.  This  includes  the  simulation  of  internal  and  external  SRM  exhaust  flows  at  lower  altitudes,  where  the  gas 
phase  may  be  accurately  modeled  using  continuum  CFD  methods.  To  the  authors’  knowledge,  this  is  the  first 
implementation  of  a  Monte  Carlo  ray  trace  model  to  allow  for  strong  coupling  between  radiation  and  flowfield 
characteristics  in  a  high  speed  emitting,  absorbing  and  scattering  multiphase  medium  of  arbitrary  optical  thickness. 
With  sufficient  modifications,  we  expect  that  the  general  procedure  described  here  may  be  applied  to  a  variety  of 
gas-particle  flows. 
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Table  1.  Particle  properties  at  the  nozzle  exit. 


Diameter, 

(im 

Mass  flux, 
kg/m2s 

Temperature, 

K 

Speed, 

m/s 

Liquid  mass 
fraction 

0.3 

0.0443 

1562 

2992 

0.579 

0.4 

0.0367 

1634 

3051 

0.661 

0.6 

0.133 

1834 

3023 

0.89 
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0.592 

2293 

2973 

1 

2 

2.29 

1920 

2855 

0.989 

4 

7.53 

2178 

2674 

1 

6 

1.84 

2407 

2472 
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Figure  1.  Contours  of  mass  density  for  particles 
and  gas.  Values  are  given  in  kg/m3. 


Figure  2.  Mean  temperature  contours  for  0.4  pm 
and  4  pm  diameter  particles.  Values  are  in  K. 
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Figure  3.  Profiles  of  particle  temperature  100  m 
downstream  of  the  nozzle  exit. 


Figure  4.  Mean  convective  and  radiative  heat 
transfer  rates  along  the  centerline. 
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Figure  5.  Contours  of  the  net  direction-averaged  Figure  6.  Comparison  of  the  net  radiative  energy 

radiative  energy  flux.  Values  are  in  W/m2.  flux  with  and  without  searchlight  emission.  Values 

are  in  W/m2. 
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Figure  7.  Energy  flux  contours  for  simulations 
employing  anisotropic  (top)  and  isotropic 
scattering  models.  Values  are  in  W/m2. 


Figure  8.  Spectral  radiance  at  a  sensor  20  cm 
upstream  of  the  nozzle  exit. 
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